## ssa -- 2.1 Stochastic Simulation Algorithm

ssa <- function(tfinal, x0, transitions)
{
	## step 1 -- initialize x0
	x <- x0
	t <- 0
	df <- data.frame()
	while(t < tfinal) {
		## step 2 -- calculate transition rates
		a <- transitions(x)
		## step 3 -- sum of transitions
		asum <- sum(a)
		## step 4 -- rexp with mean 1/asum
		tau <- rexp(1, asum)
		## step 5 -- draw from discrete distribution
		i <- which(runif(1) < cumsum(a)/asum)[1]
		## step 6 -- update time and system state
		t <- t + tau
		x <- x + V[,i]
		## step 7 -- rinse and repeat
		df <- rbind(df, data.frame(time=t, U=x[1], C=x[2], J=x[3]))
	}
	df
}
